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We investigated a nonlinear advection-diffusion-reaction equation for a passive scalar field. The 
purpose is to understand how the compressibility can affect the front dynamics and the bulk burning 
rate. We study two classes of flows: periodic shear flow and cellular flow both in the case of fast 
advection regime, analysing the system at varying the extent of compressibility and the reaction 
rate. We find that the bulk burning rate Vf in a shear flow increases with compressibility intensity, 
e, following a relation Avf ~ e 2 . Furthermore, the faster the reaction the more important the 
difference with respect to the laminar case. The effect has been quantitatively measured and it 
turns out to be generally little. For the cellular flow, the two extreme cases have been investigated, 
with the whole perturbation situated either in the centre of the vortex or in the periphery. The 
dependence in this case does not show a monotonic scaling with different behaviour in the two cases. 
The enhancing remains modest and always less than 20%. 
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I. INTRODUCTION 



The dynamics of reacting species presents several issues of great interest from a theoretical point of view [1-3]. 
Moreover it is also a problem of wide application in many fields, from front propagation in gases [4], chemical reaction 
qv • in liquids [5, 6] and ecological dynamics of biological systems (e.g. plankton in oceans) [7-11]. 

In the most simplest model of reaction dynamics, the state of the system is described by a single scalar field 
8(r,t), that represents the concentration of products. The field 8 vanishes in the regions filled with fresh material 
(the unstable phase), equals unity where only inert products are left (the stable phase) and takes intermediate 
rf\ . values wherever reactants and products coexist, i.e., in the region where production takes place. In their seminal 
contributions, Fisher, Kolmogorov, Petrovskii and Piskunov [12, 13] (FKPP) considered the simplest case of pure 
reaction/diffusion and proposed the so-called FKPP model 

><: d t 6 = DA6 + f(e), (LI) 

where D is the molecular diffusivity and f(9) describes the reaction process that obviously depends on the phenomenon 
under investigation. In this work, as in the original works of FKPP, we focus on pulled reaction, e.g. the autocatalytic 
reaction / = a0{l — 9), where a is the reaction rate and its inverse, r = 1/a, is the reaction time. 

However, most natural phenomena take place in deformable media like fluids and therefore transport properties 
cannot be ignored. If the medium is stirred, e.g. an eulerian velocity field u(x) is present, Eq. (1.1) can be generalized 
in the incompressible case to 

d t + (u • V)0 = DAO + f(0) . (1.2) 

The complete mathematical description of these phenomena is given by partial differential equations (PDE) for the 
coupled evolution of the velocity field and of the concentration of the reacting species [4]. Therefore the above 
Eq. (1.2) should be coupled with Navier-Stokes equations (usually in a non trivial way). This is the general framework 
for treating engineering combustion problems in gases [14-16]. In some cases, e.g. [17], the coupling can be simplified 
using a Boussinesq term. 



In this work, as a further simplification, we assume that the reactants do not influence the velocity field which evolves 
independently. In such a limit the dynamics is still non trivial and it is completely described by the above ARD 
Eq. (1.2) together with the proper definition of a given velocity field, u(x). This equation has been intensively studied 
in incompressible media [18-21]. In particular, it has been investigated the dependence of the front speed as a function 
of D, a and the velocity field u(x) [22]. 

On the contrary, in the case of compressible flows, the ARD problem did not receive too much attention but only 
recently in a mathematical framework [23, 24], To account for compressible flows is indeed not simple but it is 
a relevant issue in combustion [14, 16], plankton dynamics in turbulent flows [25] and also in particle-laden flows, 
where the particle phase can be highly compressible even in incompressible flows, because of inertia [26-28]. While the 
passive scalar approximation for reactive species is hardly tenable in gas combustion phenomena, it may be considered 
appropriate in aqueous or liquid reactions (notably plankton in oceans) and for dilute particle- laden flows. In those 
cases, it may give some relevant insights for front propagation and can be used as a model for the flame tracking in 
some limits [2 9]. 

Our aim is to investigate the effects of the compressibility to the bulk burning rate of the reaction process by 
studying the following PDE: 
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The scalar field 8 represents the mass fraction of a single species of a binary mixture, while Ui is the i th component 
of a given compressible flow, Do = pD is the diffusion coefficient (supposed to be constant), f(9) — u> the rate of 
production of the chosen species and p the non constant density of the fluid. 

The paper is organized as follows: Section II is devoted to the presentation of the model and the principal aspects 
of the numerical computations. In Section III we discuss the results for the front propagation in compressible shear 
flows. Section IV is devoted to the case of compressible cellular flows. Finally, in Section V, the reader can find the 
conclusions. 

II. THE MODEL 

The PDE model described by Eq. (1.3) can be derived from the equation of conservation of species, that is relevant 
for combustion dynamics [4, 15]. Let us consider two species (namely A,B) which diffuse and react together while 
they are passively transported by a compressible flow, being p A {x,y,t) the mass of species A per unit volume, the 
conservation of species A gives: 

^± + l-[ PA (u i + U A , i )}=uj A (II.l) 

where Ui is the i th component of the advective flow field, Ua,i is the velocity of diffusion of species A and u>a is the 
rate of production. 

Define the mass fraction Yk = pk/ p, where p is the density of the mixture and k = A, B. The species conservation 
can be written in terms of mass fraction as follows: 

^~ + -^ IpYk («i + U k ,i)] = u; k (II.2) 

Where Ya + Yb = 1. Moreover, if Fick's law is considered, the diffusion velocities can be defined as follow: 

8Ya 
Y A U A ,i = -Y B U B ,i = - D -Q^ (H.3) 

In the following, we assume an auto-catalytic irreversible law A + B — > 2A: 

lu a = ap A pB = ap 2 Y A Y B = ap 2 Y A (l - Y A ) (II.4) 

where the constant a controls the speed of reaction and by definition lo a = —Cob- 

Thus, the evolution of the mass fraction of species A is completely described by the following PDE: 



dt dxi 
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Dog-2+ap A 9(l-6), (II.5) 



that holds if we neglect the coupling between conservation of species equation and the conservation of energy equation. 
That is the case in which the energy released by the reaction is negligible and thus the momentum and energy equations 
evolve independently. The left hand side of Eq. (II. 5) is written in non-conservative form using the continuity equation 
of the mixture and the product pD = Do is assumed constant (which is quite a reasonable hypothesis [14, 15]). 
Since we are interested, in the front propagation, we consider the following geometry: 



-oo < x < oo , < y < L 



(11.6) 



Only for the sake of simplicity we assume periodic boundary conditions in the y— direction and 9(—oo,y,t) = 1 
(burned material in a combustion terminology) and 6(oo,y,t) = (fresh material). 
At t = the initial condition is given by: 



9(x 7 y,t) 



1 if x < 
if x > 



(II.7) 



Of course different boundary and initial conditions may be interesting. For instance, if one is interested in quenching 
issues, appropriate initial conditions would pose 6 initially localized in a region of size I 



(x,y,0) = 



1 if -1/2 <x<£/2 

if x < -1/2 on> £/2 
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Equation (II. 5) has been solved using a eighth-order central finite difference scheme in space and a fourth-order 
Runge-Kutta integration in time. The grid size is sufficiently small to guarantee a good representation of the shear 
across the reacting region and convergence of solutions has been verified. To compute accurately the asymptotic mean 
bulk burning rate, very long integration periods are required. The grid is remapped following the reacting front and 
the computational domain is extended upstream and downstream from the reactive zone so that the boundary effects 
are negligible. 



III. COMPRESSIBLE SHEAR FLOW 
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Figure 1: Shape of the active part of the front (we use the function 46(1 — 0), which is maximal for 9 = 0.5), for a fixed Peclet 

number (Pe — 100) and for two different reaction rates. Taking as reference the incompressible test case (e — 0) we can define 

two different thicknesses. The bare front thickness as 8 ~ \J D/ol and the distance between the tip and the tail of the reacting 

region A. For a slow reaction (upper panel, Da = 1), we have approximately A ~ 10 and 8 ~ 6. For fast reaction (lower 

panel, Da = 100) A « 0.8 and 8 « 0.15. 



We investigate the effects of compressibility in a 2D steady-state shear flow using the velocity field 



u(x,y) 
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(III.1) 



Such a choice corresponds to a Kolmogorov flow with amplitude Uq and wavelength L, perturbed by a steady wave 
of wavelength A and magnitude e accounting for the compressibility of the flow. Let us stress that the perturbation 
is oriented along the direction of propagation of the reactive front, i.e. the x axis. 
In order to satisfy the continuity equation, di(pui) = 0, it is necessary to impose a spatially dependence on p, as 



P(x) = Po 



1 + esm^— 



Finally, Eq. (II. 5) can be written in non-dimensional form: 

do ae i 



dt* l dx* p*Pe d. 
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p*Da6(l - 6) (111.2) 



1 <- T 



if we define p* = p/ po, x* = Xi/L, u* = Ui/Uo, t* = (tUo)/L. 

The adimcnsional parameters Pe = (poUoL) / Do and Da = (Lapo)/Uo are the Peclet and the Damkohlcr numbers 
which define the ratio between the diffusive and advective time scale and the ratio between the advective and reactive 
time scale respectively. 

In the following, we will drop the star notation and we will solve Eq. (III. 2) focusing on regimes at high Peclet 
number Pe 3> 1. Varying the Damkohler number in a range of Da G [1, 1000], we will quantify the effects of A and e 
on the asymptotic value of the bulk burning rate. 

The instantaneous bulk burning rate is: 

[•1 p-\-oo pi /> + oo 

v f (t)= / Lodxdy=l I Dap 2 6(l - 6) dxdy , (III.3) 

Jo Jo Jo Jo 

while the mean or asymptotic bulk burning rate is defined as the time average of Vf (t) over a sufficiently long interval: 

T , v f (t)dt (III A) 

To shed some light on the role played by A we first run a simulation in absence of compressibility for two different 
Damkohlcr (slow and fast reaction) and for a fixed Peclet. We characterize the thicknesses of the reactive front A 
and S (sec Figure 1 for definition). From this figure, it is clear that the faster the reaction the thinner the flame. 

Then we have carried out simulations in which the compressibility is fixed (e = 0.5) and we choose A approximately 
greater, lower or between the two thicknesses computed in the case of zero compressibility. 

In Figure 2 we show how the geometrical aspect of the reactive front changes by varying A. In the low density 
zones, the front thickness appears broader than in high density zones due to the decreasing of the local Peclet and 
Damkohlcr number (Pe; = pPe, Dai — pDa). Compressibility perturbation wrinkles the front in the small- wavelength 
limit, whereas for large wavelengths it is only corrugated, since the entire reactive-diffusive front is embedded in a 
wavelength. Nevertheless, even though the front does not appear stationary (even in a co-moving reference system) 
and it is noticeably distorted by the presence of compressibility, the mean velocity of propagation (v) does not change, 
see Figure 3. The wavelength of the perturbation controls the frequency and the magnitude of the instantaneous 
value of the front speed but does not affect the mean value. Since asymptotic propagation is not affected by A, from 
now on in all simulations we set A = 1. 

In order to quantitatively characterize the effect of compressibility, we vary both the parameter e and Da, with a 
fixed Peclet number. For this purpose, it is convenient to define the percentage difference of the mean asymptotic 
front speed between the compressible and the incompressible case as follow: 

Av % = 100 ^-^- (111.5) 

v 

where v° is the asymptotic bulk burning rate, as defined in (III. 4), for the incompressible case (e = 0). Results are 
shown in Figure 4. 

In general, in the regimes investigated here, we observe that the presence of compressibility can slightly improve the 
process of reaction and the effects grow by increasing both e and Da. For a fixed characteristic reaction rate (see 
Figure 4. a), numerical simulations suggest a power (quadratic) law of the velocity enhancement as a function of the 
parameter e 

Ati % ~£ 2 . (III.6) 
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Figure 2: Snapshot of 4(9(1 — 6) for a fixed Peclet number (Pe — 100) and for Damkohler number Da = 100. Panel (a) refers 

to a incompressible simulation (e = 0) while for the others e = 0.5. For the compressible tests the characteristic length A is set 

to be approximately lower (panel (b), A = 0.1), between (panel (c), A = 0.5) or greater (panel (d), A = 4) than the two 

thicknesses A and 8. 





(c) 

Figure 3: (Color online) Front speed and burnt mass (rribit) — L Vf(t) di) as function of time for a fixed Pe = 100, two 
different Damkohler (Da — 1 in panel: a,c ; Da — 100 in panel: b,d) and different A. 



Instead, the dependence on Damkohler is much slower. As shown in Figure 4.b the parameter Aw% is always positive 



an it grows following (approximately) a logarithmic law: 

Av% ~ a hi(Da) + b 



(III.7) 



where a and 6 may depend on e. Therefore even in the case of very strong compressibility (e = 0.5) and very fast 
reaction (a = 1000) the difference never exceeds a modest 6%. The effect of compressible wave perturbations appears 
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Figure 4: (Color online) Comparison of the front speed between compressible and incompressible shear flow for a fixed Peclet 
number (Pe = 100) at different compressibility magnitude e (on the left) and for different Damkohler (on the right). 

therefore to be: i) wrinkling; ii) second-order enhancement. 

IV. COMPRESSIBLE CELLULAR FLOW 

We discuss now the case of cellular flows, i.e., 2D steady flows of amplitude Uq composed by counter-rotating vortex 
of dimension L/2. The compressibility is imposed in the following way: 



P(x,y) = p Q C(x,y) 

'C/ oSm (HM) cos (^) _ [/oCOS (2 F)sin (2 F ) 
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(IV.l) 
(IV.2) 



We choose two different shapes for C(x,y). In the first, that we call (a) case, the density of the mixture is higher in 
the centre of the vortex: 



C(x,y) 



2ttx\ {2iry\ 
sm [ sin 



In the second, that we call (b) case, the density is higher in the periphery of the vortex: 
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(IV.3) 



(IV.4) 



The two different configurations are shown in Figure 5. The constant factor — j has been introduced in order to 
have a density perturbation which is zero in average. As in the case of the shear flow, we study the dependence of 
the dynamics on the compressibility intensity e and Damkohler number in the more realistic case of fixed high Peclet 
number. 

We will consider a wider range of Damkohler exploring the regimes at Da -c 1, Da « 1 and Da ;§> 1. Nevertheless 
we will remain in regimes PeDa > 1 which means that the characteristic diffusion time is always larger than the 
reaction time. Unlike the shear flow, in the cellular flow Av% does not show a monotonic dependence neither on e 
nor on Da, as it can be seen from Figure 6. 

Such a feature has been observed also for other configurations of C(x, y) (simulations not shown here) confirming 
that the non monotonic behaviour of Av% is related to the whirling geometry of the flow rather than to the choice 
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Figure 5: Compressible Cellular flow : (a) the density is higher in the centre of the vortex, (b) The density is higher in the 

outer region. 

of the density perturbation. In the slow reaction regime (Da <C 1) it can be seen in Fig. 6 an almost Damkohler- 
indipendent behaviour of Av% in both cases (a) and (b). On the other hand, in a middle range of Damkohler where 
the combined effect of advection and reaction is more intriguing, the two flow configurations show opposite trends 
for Av%, the reaction is faster, when the density is higher in the centre (case (a)), whereas it is slower when the 
perturbation is at the periphery (case (b)). Such a behaviour is not surprising, since the interplay between reaction 
and diffusion in the presence of closed streamlines can lead to a non trivial behaviour also in the case of uncompressible 
flow [20] , and the presence of variations in the density of the flow can act in a very non intuitive way. 
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Figure 6: (Color online) (a,b) Percentage difference between compressible and incompressible test case of the mean 
asymptotic bulk burning rate. 1/Pe — 0.003. In panel (a) the density of the fluid is higher in the centre of the vortexes while 

in panel (b) the density is higher in the periphery. 



CONCLUSIONS 



We have studied the propagation of fronts through an advection-diffusion-reaction equation where the nonlinear 
reaction term is given by the classical FKPP source term. The advective flow is generated by an imposed field which 
is perturbed by compressible waves. The compressibility is controlled by the parameter e. Two velocity fields have 
been considered: a shear-flow and a cellular one. 

In the considered flows, the front can be strongly affected by compressibility and the compressibility field forces a 
strong localization of density, but the quantitative differences with respect to the incompressible model appear modest 
(of the order of some percent). On the basis of previous studies, we do not think that the presence of chaos (turbulent 
fluctuations) should change much the scenario [21]. 



Some comments are in order to discuss the apparent difference of the behaviour of At>% in the cases of shear 
flow and cellular flow (see Figs. 4 and 6). The stream lines in the two cases are very different: namely open and 
closed, respectively. In the shear flow, the effect of compressibility on the front propagation is only slightly modified 
with respect to the uncompressible case, since the front is mainly driven by the stream. On the other hand, closed 
streamlines trigger entangled mechanisms between reaction and diffusion, that, coupled with the compressibility 
generate highly non trivial features. An example of this complicated behaviour can be found in the non monotonic 
dependence of Av% by the Damkdhler number, or in the well apparent difference between cases (a) and (b) of the 
cellular flows here considered. 

Finally, it is interesting to note that a similar model has been recently used for the study of population dynamics 
in turbulent flows [24] : 

3C 

— +V-(uC)=D \/ 2 C + fiC(l-C) (V.f) 

at 

where the scalar C(x, t) is the concentration of a population [24], which is the equivalent of our pQ in Eq. (II. 5). 
When V • u ^ 0, clustering of the population near compression regions (V • u < 0) is observed. In those regions, 
the concentration can take values greater then one and reaction rate on Eq. (V.l) can be negative, so that the scalar 
C(x, t) is not a fractional parameter. Within this model, authors linked changes in the overall carrying capacity of 
the ecosystem (i.e. the density of biological mass of the system) to the compressibility and its effect of localisation. 
Yet, present results show that the change of the carrying capacity is not due to density localisation, but rather to the 
choice of a different reaction term, which allows negative rate in high density zones. Indeed, in the present work we 
have a strong density localisation but our FKPP model for a fractional parameter does not allow negative rate. The 
results is that the average carrying capacity does not change even in presence of compressible flows. The analysis 
of present results for the Lagrangian displacement of passive reactive tracers and irreversible reaction dynamics is 
ongoing. 
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